### This project replicates the analyses reported in Toshkov and Krouwel (2022) 'Beyond the U-Curve: Citizen Preferences on European Integration in Multidimensional Political Space' in European Union Politics 
### Part I: This script prepares and analyzes the data from the 2019 EP elections in The Netherlands
# Clean up VAA data and make variables -------------------------------------------------------------------------
## This part of the script shows the preparation of the datafile that is shared from the original datafile that is proprietary
# Load the data
nl.all <- read_sav('./data in/nl ep 2019/2019 VAA data AN.sav')

# Recode variables and filter
nl.all<- nl.all %>%
  mutate(lr = popup_answer_1040 - 5,
         proeu = popup_answer_1041 - 5,
         progr = popup_answer_1042 - 5,
         ep2014.vote = recode (popup_answer_1039, '0' = 'D66', '1' = 'CDA', '2' = 'PVV', '3' = 'VVD', '4' = 'SP', '5' = 'PvdA', 
                               '6' = 'CU', '7' = 'GL', '8' = 'PvdD', '9' = 'Other', '10' = 'Other', 
                               '11' = NA_character_, '12' = 'No vote', '13' = NA_character_),
         ep2019.vote.intent = recode (popup_answer_1038, '0' = '50+', '1' = 'CDA', '2' = 'CU', '3' = 'D66', '4' = 'Other', '5' = 'DENK', 
                               '6' = 'FvD', '7' = 'GL', '8' = 'Other', '9' = 'PvdA', '10' = 'PvdD', '11' = 'PVV', '12' = 'SP', 
                               '13' = 'Other', '14' = 'Other',  '15' = 'VVD',
                               '16' = NA_character_, '17' = 'No vote', '18' = NA_character_),
         tk2017.vote = recode (popup_answer_1036, '0' = 'VVD', '1' = 'PVV', '2' = 'CDA', '3' = 'D66', '4' = 'GL', '5' = 'SP', 
                                      '6' = 'PvdA', '7' = 'CU', '8' = 'PvdD', '9' = '50+', '10' = 'SGP', '11' = 'DENK', '12' = 'FvD', 
                                      '13' = 'Other', '14' = NA_character_,  '15' = 'No vote',
                                      '16' = NA_character_, '17' = NA_character_),
         vote2017 = recode (popup_answer_1036, '0' = 'VVD', '1' = 'PVV', '2' = 'CDA', '3' = 'D66', '4' = 'GL', '5' = 'SP', 
                               '6' = 'PvdA', '7' = 'CU', '8' = 'PvdD', '9' = '50+', '10' = 'SGP', '11' = 'DENK', '12' = 'FvD', 
                               '13' = NA_character_, '14' = NA_character_,  '15' = NA_character_,
                               '16' = NA_character_, '17' = NA_character_),
         ep.2019.vote.sure = recode (popup_answer_1037, '0' = 'Sure', '1' = 'Doubt', '2' = 'Dont_know', '3' = 'No vote'),
         age = as.numeric(as.character(ifelse (bg_question_445==94, NA_character_, 2019 - (2003 - bg_question_445)))),
         age_cat3 = ifelse(age<35, 'young', ifelse(age<65, 'midage', 'old')),
         sex = recode (bg_question_416, '0' = 'man', '1' = 'woman', '2' = NA_character_, '3' = NA_character_),
         edu = recode (bg_question_447, '0' = '4.high', '1' = '4.high', '2' = '3.midhigh', '3' = '2.midlow', '4' = '2.midlow', '5' = '1.lower',
                       '6' = '1.lower', '7' = NA_character_),
         edu_cat3 = ifelse(edu=='4.high', 'high_edu', ifelse(edu=='3.midhigh', 'mid_edu', 'low_edu')),
         themes = case_when ((relevant_themes_130 == 1 ~ 'society'),  (relevant_themes_162 == 1 ~ 'security'),
                             (relevant_themes_268 == 1 ~ 'eu_integration'), (relevant_themes_269 == 1 ~ 'immigration'),
                             (relevant_themes_271 == 1 ~ 'work_incomes'), (relevant_themes_273 == 1 ~ 'economy'),
                             (relevant_themes_274 == 1 ~ 'environment')),
         sure = recode (popup_answer_1037, '0' = 1, '1' = 0, '2' = 0, '3' = 0, '4' = 0),
         type = 'voter',
         index = paste(sex, age_cat3, edu_cat3, sep='.'),
         pindex = paste(vote2017, sex, age_cat3, edu_cat3, sep='.'),
         eu.bad = statement_answer_4907 - 3,
         privacy.limit = statement_answer_4928 - 3,
         eu.money.climate = statement_answer_4923 - 3,
         many.immigrants = statement_answer_4938 - 3,
         euro.out = statement_answer_4908 - 3,
         eu.army = statement_answer_4926 -3,
         eu.money.culture = 6 - (statement_answer_4925) -3,
         weed.legal = statement_answer_4913 -3,
         eu.enforce = statement_answer_4932 -3,
         eu.money.aid = statement_answer_4941 -3,
         greener.energy = statement_answer_4929 -3,
         eu.nat.borders = statement_answer_4921 -3,
         eu.money.sustainable = statement_answer_4927 -3,
         soc.sec.stable = statement_answer_4934 -3,
         islam.bad = statement_answer_4920 -3,
         market.health = statement_answer_4918 -3,
         homo.adopt = statement_answer_4912 -3,
         eu.veto = statement_answer_4911 -3,
         fire.easy = statement_answer_4917 -3,
         eu.fp = statement_answer_4910 -3,
         state.econ.out = statement_answer_4916 -3,
         eu.lost.id = statement_answer_4922 -3,
         eu.fin.solidarity = statement_answer_4924 -3,
         limit.asylum = statement_answer_4935 -3,
         eu.unempl = statement_answer_4937 -3,
         redistribute = statement_answer_4915 -3,
         fighters.back = statement_answer_4939 -3,
         eu.out = statement_answer_4940 -3,
         eu.tax = statement_answer_4931 -3,
         eu.limit.work = statement_answer_4909 -3
         ) %>%
  filter (country == 'Netherlands', amount_visits == 1,
          seconds_on_statements > 30*3, seconds_on_statements < 30*30) %>%
  select (lr:ncol(.))

# remove those who answered everything with Neutral
nl.all$zeroes <- nl.all %>% 
  select (eu.bad:eu.limit.work)  %>%
  abs () %>%
  rowSums(na.rm=T)

nl.all <- nl.all %>% filter (zeroes!=0) %>% select (-zeroes)

# create deductively derived positions on scales
nl.all <- nl.all %>%
  mutate (obj.lr = (- redistribute - soc.sec.stable + fire.easy + state.econ.out + market.health)/5,
          obj.cp = (- islam.bad - privacy.limit - limit.asylum - many.immigrants +  homo.adopt +
                      weed.legal + greener.energy + fighters.back)/8,
          obj.eu = (- eu.bad - euro.out - eu.veto - eu.lost.id - eu.out + eu.army + eu.tax + 
                      eu.enforce + eu.fp + eu.unempl)/10,
          obj.eu2 = (- eu.bad - euro.out - eu.veto - eu.lost.id - eu.out + eu.army + eu.tax + 
                      eu.enforce + eu.fp - eu.limit.work - eu.nat.borders + eu.money.climate +
                      eu.money.sustainable + eu.unempl + eu.fin.solidarity + eu.money.aid + eu.money.culture)/17
          )


nl.all$id = seq.int(nrow(nl.all)) #index for merging later

# save
save (nl.all, file = './data out/nl2019_replication.RData')

# Load and subset the data (start replication here) -------------------------------------------------------------------------
load (file = './data out/nl2019_replication.RData')

# Subset policy items only
nl.s <- nl.all %>% select(id, eu.bad:eu.limit.work) %>%
  filter(complete.cases(.))

# Number of factors ----------------------------------------------------------------
# get polychoric covariance
tk.poly<-polychoric(select(nl.s, -id))

# determine number of factors
ev <- eigen(cor(select(nl.s, -id)))
ev.poly <- eigen(tk.poly$rho) #with the polychronic correlations
# number of factors
ap <- parallel(subject=nrow(nl.s),var=ncol(select(nl.s, -id)),rep=10,quantile=.95)
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea)
plotnScree(nS)

# PCA analysis (Table A5) ----------------------------------------------------------------
pca.1 <- principal(select(nl.s, -id), nfactors=5, rotate='none', cor='poly', scores=TRUE)
print(pca.1, digits=2, cut=.31, sort=T)

pca.out<-as.data.frame(unclass(pca.1$loadings))
pca.out <- pca.out %>%
  mutate_at(vars(starts_with('PC')), ~replace(., abs(.) < 0.4, 0)) %>%
  mutate_at(vars(starts_with('PC')), funs(round(., 2))) %>%
  arrange (-abs(PC1),-abs(PC2), -abs(PC3), -abs(PC4), -abs(PC5)) %>%
  mutate_at(vars(starts_with('PC')), ~replace(., abs(.) == 0, '')) %>%
  mutate ( issue = rownames(.)) %>%
  relocate (issue) 
pca.out <- rbind (pca.out, c('Proportion Variance', paste0(round(prop.table(pca.1$values)[1:5]*100,0), '%')))
pca.out

pca.nl.table = htmlTable::htmlTable(pca.out, rnames = FALSE,
                                    col.columns = c('white','lightgrey'),
                                    tfoot = 'Note: Loadings smaller than +/-0.41 are not printed',
                                    align = c('l','r','r','r','r','r'),
                                    align.header = c('l','r','r','r','r','r'),
                                    header = c('Policy issue','PC1','PC2','PC3','PC4','PC5'),
                                    caption = "Table X. PCA results from the Netherlands, 2019")
print(pca.nl.table, type = "html", file = './tables/pca_nl_table.html')


# Factor analysis I (full sample) (Table 1)----------------------------------------------------------------
fa.1 <- fa (select(nl.s, -id), nfactors=6, rotate='Promax', cor='poly', scores=TRUE)
print(fa.1, digits=2, cut=.33, sort=T)

fa.out<-as.data.frame(unclass(fa.1$loadings))
fa.out <- fa.out %>%
  mutate_at(vars(starts_with('MR')), ~replace(., abs(.) < 0.31, 0)) %>%
  mutate_at(vars(starts_with('MR')), funs(round(., 2))) %>%
  arrange (-abs(MR2),-abs(MR1), -abs(MR3), -abs(MR5), -abs(MR4), -abs(MR6)) %>%
  mutate_at(vars(starts_with('MR')), ~replace(., abs(.) == 0, '')) %>%
  mutate ( issue = rownames(.)) %>%
  relocate (issue) 
fa.out <- rbind (fa.out, c('Proportion Variance', paste0(round(fa.1$Vaccounted[2,]*100,0), '%')))

fa.out

fa.nl.table = htmlTable::htmlTable(fa.out, rnames = FALSE,
                                   col.columns = c('white','lightgrey'),
                                   tfoot = 'Note: Loadings smaller than +/-0.31 are not printed',
                                   align = c('l','r','r','r','r','r','r'),
                                   align.header = c('l','r','r','r','r','r','r'),
                                   header = colnames(fa.out),
                                   caption = "Table X. Factor Analysis results from the Netherlands, 2019")
print(fa.nl.table, type = "html", file = './tables/fa_nl_table.html')

# Factor analysis II (restricted sample) (Table A6)----------------------------------------------------------------

# Run the same factor model but only on the subset who filled in the additional module
nl.s2 <- nl.all %>%
  select (id, lr, progr, proeu, obj.lr, obj.cp, obj.eu, obj.eu2, eu.bad:eu.limit.work) %>%
  na.omit 

# factors with fa.poly and promax rotation
fa.2 <- fa (select(nl.s2, eu.bad:eu.limit.work), nfactors=6, rotate='Promax', cor='poly', scores=TRUE)
print(fa.2, digits=2, cut=.33, sort=T)

fa2.out<-as.data.frame(unclass(fa.2$loadings))
fa2.out <- fa2.out %>%
  mutate_at(vars(starts_with('MR')), ~replace(., abs(.) < 0.31, 0)) %>%
  mutate_at(vars(starts_with('MR')), funs(round(., 2))) %>%
  arrange (-abs(MR4),-abs(MR1), -abs(MR3), -abs(MR5), -abs(MR2), -abs(MR6)) %>%
  mutate_at(vars(starts_with('MR')), ~replace(., abs(.) == 0, '')) %>%
  mutate ( issue = rownames(.)) %>%
  relocate (issue) 
fa2.out <- rbind (fa2.out, c('Proportion Variance', paste0(round(fa.2$Vaccounted[2,]*100,0), '%')))

fa2.out

fa2.nl.table = htmlTable::htmlTable(fa2.out, rnames = FALSE,
                                    col.columns = c('white','lightgrey'),
                                    tfoot = 'Note: Loadings smaller than +/-0.31 are not printed',
                                    align = c('l','r','r','r','r','r','r'),
                                    align.header = c('l','r','r','r','r','r','r'),
                                    header = colnames(fa2.out),
                                    caption = "Table X. Factor Analysis results from the Netherlands, 2019 [interested sample]")
print(fa2.nl.table, type = "html", file = './tables/fa2_nl_table.html')

# Factor analysis III (only those who are sure) (Table A7----------------------------------------------------------------
# Run the same factor model but only on the subset who are sure who they will vote for
nl.s4 <- nl.all %>%
  filter(sure == 1) %>%
  select (eu.bad:eu.limit.work) %>%
  na.omit 

# factors with fa.poly and promax rotation
fa.4 <- fa (select(nl.s4, eu.bad:eu.limit.work), nfactors=6, rotate='Promax', cor='poly', scores=TRUE)
print(fa.4, digits=2, cut=.33, sort=T)

fa4.out<-as.data.frame(unclass(fa.4$loadings))
fa4.out <- fa4.out %>%
  mutate_at(vars(starts_with('MR')), ~replace(., abs(.) < 0.31, 0)) %>%
  mutate_at(vars(starts_with('MR')), funs(round(., 2))) %>%
  arrange (-abs(MR1),-abs(MR5), -abs(MR3), -abs(MR6), -abs(MR2), -abs(MR4)) %>%
  mutate_at(vars(starts_with('MR')), ~replace(., abs(.) == 0, '')) %>%
  mutate ( issue = rownames(.)) %>%
  relocate (issue) 
fa4.out <- rbind (fa4.out, c('Proportion Variance', paste0(round(fa.4$Vaccounted[2,]*100,0), '%')))

fa4.out

fa4.nl.table = htmlTable::htmlTable(fa4.out, rnames = FALSE,
                                    col.columns = c('white','lightgrey'),
                                    tfoot = 'Note: Loadings smaller than +/-0.31 are not printed',
                                    align = c('l','r','r','r','r','r','r'),
                                    align.header = c('l','r','r','r','r','r','r'),
                                    header = colnames(fa4.out),
                                    caption = "Table X. Factor Analysis results from the Netherlands, 2019 [certain sample]")
print(fa4.nl.table, type = "html", file = './tables/fa4_nl_table.html')

# Assign factor loadings to cases ----------------------------------------------------------------
nl.s <- bind_cols(nl.s, data.frame(pca.1$scores), data.frame(fa.1$scores))
nl.all <- left_join (nl.all, select(nl.s, id, PC1:MR6), by='id')

# Subset -------------------------------------------------
nl.s <- nl.all %>% 
  filter (type=='voter') %>% 
  select ('lr','progr','proeu') %>% 
  filter(complete.cases(.))

# Self-placement analysis: 3D plots -------------------------------------------------
#prepare data
elevation.df = data.frame(x = nl.s$lr, y = nl.s$progr, z = nl.s$proeu)
elevation.loess = loess(z ~ x*y, data = elevation.df, degree=2, span=0.4)
xmin<- -5
xmax<- 5
ymin<- -5
ymax<- 5

elevation.fit = expand.grid(list(x = seq(xmin, xmax, 0.1), y = seq(ymin, ymax, 0.1)))
z = predict(elevation.loess, newdata = elevation.fit)
elevation.fit$Height = as.numeric(z)

#prepare matrix with x and y and rows/cols and z as the cell values
u1<-unique(elevation.fit$x)
u2<-unique(elevation.fit$y)
temp.data<-matrix(NA, nrow=length(u1), ncol=length(u2))

for (i in 1:length(u1)){
  for (j in 1:length(u2)){
    temp.data[i,j]<-elevation.fit$Height[elevation.fit$x==u1[i] & elevation.fit$y==u2[j]]
  }
}

ix <- seq(1, nrow(temp.data), length.out = 11)
iy <- seq(1, ncol(temp.data), length.out = 11)

# 3D in black and white (Figure 2, left panel)
png('./figures/nl 2019 eu3.png', width=900*2, height=700*2, res=96)
par(mfrow = c(1, 1), mar = c(0, 0, 0, 0))
persp3D(z = temp.data, theta = 45, phi = 10, facets = T, curtain = TRUE, 
        shade=0.8, col = ramp.col(c("white", "black")), border = "grey", colkey = FALSE, xlab="Left-Right", ylab="Conservative-Progressive", zlab="Anti-Pro EU")
dev.off()

# 2D plot with the plot3d package (Figure A1 top)
par(mfrow = c(1, 1), mar = c(6, 8, 6, 4))
#first in 2D. quite nice actually
png('./figures/nl 2019 eu.png', width=900*2, height=700*2, res=96)
image2D(temp.data, clab = "Anti-Pro EU",xlab="Economic Left/Right", ylab="Conservative/Progressive",
        cex.lab=1.5, cex.axis=1, cex.main=1.5, cex.sub=2, pt.cex = 2)
dev.off()

# 4 3D plots in one figure (Figure A1 bottom)
png('./figures/nl 2019 eu6.png', width=900*2, height=700*2, res=96)
par(mfrow = c(2, 2), mar = c(3, 3, 3, 3))

persp3D(z = temp.data, theta = 55, phi = 20, facets = T, curtain = TRUE, 
        shade=0.8, colkey = FALSE, xlab="Left-Right", ylab="Conservative-Progressive", zlab="Anti-Pro EU", box=T, bty='b2',ticktype = "detailed")

persp3D(z = temp.data, theta = 45, phi = 10, facets = T, curtain = TRUE, 
        shade=0.8, col = ramp.col(c("white", "black")), border = "grey", colkey = FALSE, xlab="Left-Right", ylab="Conservative-Progressive", zlab="Anti-Pro EU")

ribbon3D(z = temp.data[ix, ], along = "y", cex.axis = 0.8, 
         curtain = T, space = 0.7, theta = 25, phi = 5, shade=0.2, colkey = FALSE, xlab="Left-Right", ylab="Conservative-Progressive", zlab="Anti-Pro EU", 
         box=T, bty='b2', ticktype = "detailed")

hist3D(z = temp.data[ix,iy], theta = 55, phi = 20, shade=0.2, colkey = FALSE, box=T, bty='b2',ticktype = "detailed",
       xlab="Left-Right", ylab="Conservative-Progressive", zlab="Anti-Pro EU")

dev.off()

# with the plotly library
plot_ly(x = u1, y=u2, z=temp.data) %>%
  layout(
    title = "Citizen positions on political scales, the Netherlands, 2019",
    scene = list(
      xaxis = list(title = "Economic Left-Right"),
      yaxis = list(title = "Conservative-Progressive"),
      zaxis = list(title = "Anti-Pro EU")
    )) %>% add_surface(contours = list(
      z = list(
        show=TRUE,
        usecolormap=TRUE,
        highlightcolor="#ff0000",
        project=list(z=TRUE)
      )
    ))

# now with the rayshader package (Figure A2)
temp.data2 = temp.data*300

raymat = ray_shade(temp.data2)
ambmat = ambient_shade(temp.data2)

png('./figures/raysahder1.png', width=1000, height=800, res=96)

temp.data2 %>%
  sphere_shade(texture = "desert") %>%
  add_water(detect_water(temp.data2), color="desert") %>%
  add_shadow(ray_shade(temp.data2,zscale=3,maxsearch = 300),0.5) %>%
  add_shadow(ambmat,0.5) %>%
  plot_3d(temp.data2,zscale=35,fov=0,theta=-197,zoom=0.75,phi=10, windowsize = c(1000,800),
          water=TRUE, waterdepth = 0, wateralpha = 0.5,watercolor = "lightblue",
          waterlinecolor = "white",waterlinealpha = 0.5)


render_label(temp.data2, x=1,y=1, z=1100, zscale=35,relativez=FALSE, text = "Left, Conservative",textsize = 1.7,linewidth = 2,freetype = FALSE,
             linecolor = 'darkred', textcolor = "darkred") 
render_label(temp.data2, x=101,y=1, z=1200, zscale=35,relativez=FALSE,
             linecolor = 'darkblue', textcolor = "darkblue", text = "Right, Conservative",textsize = 1.7,linewidth = 2,freetype = FALSE)
render_label(temp.data2, x=1,y=101, z=1400, zscale=35,relativez=FALSE,
             linecolor = 'darkred', textcolor = "darkred", text = "Left, Progressive",textsize = 1.7,linewidth = 2,freetype = FALSE)
render_label(temp.data2, x=101,y=101, z=1200, zscale=35, relativez=FALSE,
             linecolor = 'darkblue',textcolor = "darkblue", text = "Right, Progresive",textsize = 1.7,linewidth = 2,freetype = FALSE)

render_label(temp.data2, x=26,y=101, z=1400, zscale=35, relativez=FALSE, dashed= TRUE,
             linecolor = 'black', textcolor = "black", text = "Peak of EU support",textsize = 1.5,linewidth = 2,freetype = FALSE, offset=0)
render_label(temp.data2, x=85,y=15, z=300, zscale=35, relativez=FALSE, dashed= TRUE,
             linecolor = 'black', textcolor = "black", text = "Lake of Eurosceptics",textsize = 1.5, linewidth = 2,freetype = FALSE)

render_snapshot()
dev.off()

# Additional 2D and 3D density plots --------------------------------------------------
# a combination of 2d plots (Figure A3, top)
p1 <- ggplot(nl.s, aes(x=lr, y=progr) ) +
  stat_density_2d(aes(fill = ..level..), geom = "polygon", colour="white") +
  theme_bw() + xlim(-5,5) + ylim(-5,5) + xlab('Left-Right') + ylab('Conservative-Progressive') + theme(legend.position="none")

p2 <- ggplot(nl.s, aes(x=lr, y=proeu) ) +
  stat_density_2d(aes(fill = ..level..), geom = "polygon", colour="white") +
  theme_bw() + xlim(-5,5) + ylim(-5,5) + xlab('Left-Right') + ylab('Anti-Pro EU') + theme(legend.position="none")

p3 <- ggplot(nl.s, aes(x=progr, y=proeu) ) +
  stat_density_2d(aes(fill = ..level..), geom = "polygon", colour="white") +
  theme_bw() + xlim(-5,5) + ylim(-5,5) + ylab('Anti-Pro EU') + xlab('Conservative-Progressive') + theme(legend.position="none")

(p1 + plot_spacer())/(p2 +p3)

# 3D histogram with the density mapped as color (Figure A3, bottom)
smooth <- MASS::kde2d(nl.s$lr, nl.s$progr, lims = c(-5, 5, -5, 5), n = 11)$z 

par(mfrow = c(2, 2), mar = c(3, 3, 3, 3))
color = rev(rainbow(10, start = 0/6, end = 4/6))

hist3D(z = temp.data[ix,iy], theta = 35, phi = 20, shade=0.2, colvar = smooth, col=color, colkey = FALSE,
       box=T, bty='b2',ticktype = "detailed",xlab="Left-Right", ylab="Conservative-Progressive", zlab="Anti-Pro EU")

hist3D(z = temp.data[ix,iy], theta = 35 + 180, phi = 20, shade=0.2, colvar = smooth, col=color,
       box=T, bty='b2',ticktype = "detailed",xlab="Left-Right", ylab="Conservative-Progressive", zlab="Anti-Pro EU")

hist3D(colvar = temp.data[ix,iy], theta = 35, phi = 20, shade=0.2, z = smooth, col=color, colkey = FALSE,
       box=T, bty='b2',ticktype = "detailed",xlab="Left-Right", ylab="Conservative-Progressive", zlab="Density")

hist3D(colvar = temp.data[ix,iy], theta = 35 + 180, phi = 20, shade=0.2, z = smooth, col=color,
       box=T, bty='b2',ticktype = "detailed",xlab="Left-Right", ylab="Conservative-Progressive", zlab="Density")


# Subjective placement and objective positions: correlation plots (Figure 1) ----------------------------------
ggpairs(
  nl.all [, c('obj.lr', 'obj.cp', 'obj.eu', 'lr','progr','proeu')],
  upper = list(continuous = "cor"),
  #upper = list(continuous = "density"),
  lower = list(continuous = "trends"),
  diag  = list(continuous = "barDiag"),
  axisLabels = 'show',
  columnLabels = c('Obj. Left-Right','Obj. Cons.-Progr.','Obj. Anti-Pro EU',
                   'Left-Right','Conservative-Progressive','Anti-Pro EU')) + 
  theme_bw()

# Objective positions: 3D plots ----------------------------------
nl.s2 <- nl.all %>% 
  filter (type=='voter') %>% 
  select ('obj.lr','obj.cp','obj.eu') %>% 
  filter(complete.cases(.))

#prepare data
elevation.df = data.frame(x = nl.s2$obj.lr, y = nl.s2$obj.cp, z = nl.s2$obj.eu)
elevation.loess = loess(z ~ x*y, data = elevation.df)

xmin<- -2
xmax<- 2
ymin<- -2
ymax<- 2

elevation.fit = expand.grid(list(x = seq(xmin, xmax, 0.1), y = seq(ymin, ymax, 0.1)))
z = predict(elevation.loess, newdata = elevation.fit)
elevation.fit$Height = as.numeric(z)

#prepare matrix with x and y and rows/cols and z as the cell values
u1<-unique(elevation.fit$x)
u2<-unique(elevation.fit$y)
temp.data<-matrix(NA, nrow=length(u1), ncol=length(u2))

for (i in 1:length(u1)){
  for (j in 1:length(u2)){
    temp.data[i,j]<-elevation.fit$Height[elevation.fit$x==u1[i] & elevation.fit$y==u2[j]]
  }
}

ix <- seq(1, nrow(temp.data), length.out = 11)
iy <- seq(1, ncol(temp.data), length.out = 11)

smooth <- MASS::kde2d(nl.s2$obj.lr, nl.s2$obj.cp, lims = c(-2, 2, -2, 2), n = 11)$z 

# 4 3D plots (Figure A4) 
par(mfrow = c(2, 2), mar = c(1, 1, 1, 1))
color = rev(rainbow(10, start = 0/6, end = 4/6))

hist3D(z = temp.data[ix,iy], theta = 35, phi = 20, shade=0.2, colvar = smooth, col=color, colkey = FALSE, clab='Density',
       box=T, bty='b2',ticktype = "detailed",xlab="Left-Right", ylab="Conservative-Progressive", zlab="Anti-Pro EU")

hist3D(z = temp.data[ix,iy], theta = 35 - 90, phi = 20, shade=0.2, colvar = smooth, col=color, clab='Density',
       box=T, bty='b2',ticktype = "detailed",xlab="Left-Right", ylab="Conservative-Progressive", zlab="Anti-Pro EU")

hist3D(colvar = temp.data[ix,iy], theta = 35, phi = 20, shade=0.2, z = smooth, col=color, colkey = FALSE, clab='Anti-Pro EU',
       box=T, bty='b2',ticktype = "detailed",xlab="Left-Right", ylab="Conservative-Progressive", zlab="Density")

hist3D(colvar = temp.data[ix,iy], theta = 35 - 90, phi = 20, shade=0.2, z = smooth, col=color, clab='Anti-Pro EU',
       box=T, bty='b2',ticktype = "detailed",xlab="Left-Right", ylab="Conservative-Progressive", zlab="Density")

# 3D plots in black and white (Figure 2, right panel)
png('./figures/nl 2019 obj eu3.png', width=900*2, height=700*2, res=96)
par(mfrow = c(1, 1), mar = c(0, 0, 0, 0))
persp3D(z = temp.data, theta = 45, phi = 10, facets = T, curtain = TRUE, 
        shade=0.8, col = ramp.col(c("white", "black")), border = "grey", colkey = FALSE, xlab="Obj. Left-Right", ylab="Obj. Conservative-Progressive", zlab="Obj. Anti-Pro EU")
dev.off()


# Preparing the LISS data ----------------------------------

# These files are downloaded from the LISS data archive. The files are public but access requires registration
liss_pol <- read_sav('./data in/nl ep 2019/LISS/cv20l_EN_1.0p.sav')
liss_gen <- read_sav('./data in/nl ep 2019/LISS/avars_201912_EN_1.0p.sav')

liss <- left_join (liss_pol, liss_gen, by='nomem_encr')

liss <- liss %>%
  filter (leeftijd > 17) %>%
  mutate (age = leeftijd,
          age_cat = lftdcat,
          age_cat3 = ifelse(age<35, 'young', ifelse(age<65, 'midage', 'old')),
          edu = dplyr::recode (as.character(oplzon), '1'='1.lower','2'='2.midlow','3'='2.midlow','4'='3.midhigh', '5'='3.midhigh', '6'='4.high', .default=NA_character_),
          edu_cat3 = ifelse(edu=='4.high', 'high_edu', ifelse(edu=='3.midhigh', 'mid_edu', 'low_edu')),
          sex = ifelse(geslacht==2, 'woman', 'man'),
          vote2017 = dplyr::recode(as.character(cv20l307), '1' = 'VVD', '2' ='PVV', '3'='CDA', '4'='D66', '5'='GL', '6'='SP','7'='PvdA', '8'='CU', '9'='PvdD', '10'='50+',
                                   '11'='SGP', '12'='DENK', '13'='FvD',.default = NA_character_),
          index = paste(sex, age_cat3, edu_cat3, sep='.'),
          pindex = paste(vote2017, sex, age_cat3, edu_cat3, sep='.')
  ) %>%
  select (age, age_cat, age_cat3, edu, edu_cat3, sex, vote2017, index) %>%
  filter (is.na(edu_cat3)==F)

# save the weigths from LISS
save (liss, file = './data out/liss.RData')

# Factor analysis V (weighted based on demographics and party)----------------------------------------------------------------
# load the LISS weights
load (file = './data out/liss.RData')

# subset the VAA data
nl.nona <- nl.all %>%
  select (index, pindex, vote2017, sex, age_cat3, edu_cat3, eu.bad:eu.limit.work) %>%
  na.omit

# get the relative weights
w<-round(prop.table(table(liss$index, liss$vote2017), 2), 2)/
  round(prop.table(table(nl.nona$index, nl.nona$vote2017), 2), 2)

w <- melt(data.frame(w)) %>%
  mutate (pindex = paste(Var2, Var1, sep='.')) %>%
  mutate_if(is.numeric, list(~na_if(., Inf))) %>%
  select (pindex, value)

# attach back to data
nl.nona<- left_join (nl.nona, w, by='pindex')

# Factor analysis IV (weighted based on party choice and demographics) (Table A8)----------------------------------------------------------------
nl.s5 <- nl.nona %>%
  select (value, eu.bad:eu.limit.work) %>%
  na.omit 

# factors with fa.poly and promax rotation
fa.5 <- fa (select(nl.s5, -value), nfactors=6, rotate='Promax', cor='poly', scores=TRUE, weight=nl.s5$value)
print(fa.5, digits=2, cut=.33, sort=T)

fa5.out<-as.data.frame(unclass(fa.5$loadings))
fa5.out <- fa5.out %>%
  mutate_at(vars(starts_with('MR')), ~replace(., abs(.) < 0.31, 0)) %>%
  mutate_at(vars(starts_with('MR')), funs(round(., 2))) %>%
  arrange (-abs(MR4),-abs(MR1), -abs(MR3), -abs(MR5), -abs(MR2), -abs(MR6)) %>%
  mutate_at(vars(starts_with('MR')), ~replace(., abs(.) == 0, '')) %>%
  mutate ( issue = rownames(.)) %>%
  relocate (issue) 
fa5.out <- rbind (fa5.out, c('Proportion Variance', paste0(round(fa.5$Vaccounted[2,]*100,0), '%')))

fa5.out

fa5.nl.table = htmlTable::htmlTable(fa5.out, rnames = FALSE,
                                    col.columns = c('white','lightgrey'),
                                    tfoot = 'Note: Loadings smaller than +/-0.31 are not printed',
                                    align = c('l','r','r','r','r','r','r'),
                                    align.header = c('l','r','r','r','r','r','r'),
                                    header = colnames(fa4.out),
                                    caption = "Table X. Factor Analysis results from the Netherlands, 2019 [weighted demo+party]")
print(fa5.nl.table, type = "html", file = './tables/fa5_nl_table.html')

# Factor analysis V (weighted based on demographics only) (Table A9)----------------------------------------------------------------
nl.nona <- nl.all %>%
  select (index, pindex, vote2017, sex, age_cat3, edu_cat3, eu.bad:eu.limit.work) %>%
  na.omit

w<-round(prop.table(table(liss$index)), 2) /
  round(prop.table(table(nl.nona$index)), 2)

w <- melt(data.frame(w)) %>%
  mutate (index = paste(Var1, sep='')) %>%
  mutate_if(is.numeric, list(~na_if(., Inf))) %>%
  select (index, value)

# attach back to data
nl.nona<- left_join (nl.nona, w, by='index')

# Run the same factor model but only on the subset who are sure who they will vote for
nl.s6 <- nl.nona %>%
  select (value, eu.bad:eu.limit.work) %>%
  na.omit 

# factors with fa.poly and promax rotation
fa.6 <- fa (select(nl.s6, -value), nfactors=6, rotate='Promax', cor='poly', scores=TRUE, weight=nl.s6$value)
print(fa.6, digits=2, cut=.33, sort=T)

fa6.out<-as.data.frame(unclass(fa.6$loadings))
fa6.out <- fa6.out %>%
  mutate_at(vars(starts_with('MR')), ~replace(., abs(.) < 0.31, 0)) %>%
  mutate_at(vars(starts_with('MR')), funs(round(., 2))) %>%
  arrange (-abs(MR4),-abs(MR1), -abs(MR3), -abs(MR5), -abs(MR2), -abs(MR6)) %>%
  mutate_at(vars(starts_with('MR')), ~replace(., abs(.) == 0, '')) %>%
  mutate ( issue = rownames(.)) %>%
  relocate (issue) 
fa6.out <- rbind (fa6.out, c('Proportion Variance', paste0(round(fa.6$Vaccounted[2,]*100,0), '%')))

fa6.out

fa6.nl.table = htmlTable::htmlTable(fa6.out, rnames = FALSE,
                                    col.columns = c('white','lightgrey'),
                                    tfoot = 'Note: Loadings smaller than +/-0.31 are not printed',
                                    align = c('l','r','r','r','r','r','r'),
                                    align.header = c('l','r','r','r','r','r','r'),
                                    header = colnames(fa6.out),
                                    caption = "Table X. Factor Analysis results from the Netherlands, 2019 [weighted demo only]")
print(fa6.nl.table, type = "html", file = './tables/fa6_nl_table.html')

### THE END
